#include "HY1C_out.h"
/* ========================================================================== */
/* module polcor.c - functions to read and apply polarization correction      */
/*                                                                            */
/* Written By: B. Franz, SAIC, November 2005.                                 */
/*     W. Robinson, SAIC, 29 May 2009  generalize for use with VIIRS          */
/* ========================================================================== */

#include "l12_proto.h"
#include "hdf_utils.h"

#define NMIR  2  /* always 2 mirror sides                    */
#define NDET  40 /* big enough for largest detector array    */
#define NANG  7  /* big enough for largest set of table AOI or scan angles */

void polcor(l1str *l1rec, int32_t ip)
{
    static int firstCall = 1;
    static float  m12[NBANDS][NMIR][NANG][NDET];
    static float  m13[NBANDS][NMIR][NANG][NDET];
    static float  detfac[NBANDS];
    static float  maxpix;
    static int32  nang;
    static float *ang  = NULL, margin_s;
    static int32_t   nir_s;
    static int32_t   nir_l;
    
    static double *dm12;
    static double *dm13;
    static int have_xcal[NBANDS];

    static float *pxwt;  /* weights for each pixel in the line and */
    static char *pxangix;  /* index of low angle storage */

    float *polcor = l1rec->polcor;
    float *dpol   = l1rec->dpol;

    int32_t  pixnum = l1rec->pixnum[ip];
    int32_t  detnum = l1rec->detnum;
    int32_t  mside  = l1rec->mside;
//    int32_t  nbands = l1rec->nbands;
    int32_t  nbands = 8;
    int32_t  sensorID = l1rec->sensorID;

    float m1 [NANG];
    float alpha;
    int32_t  ib, ipb;
    int32_t  idet, iang;
    float wt, zang;
    float L_x, L_qp, L_up;
    
    int iagsm[] = {0,640,1376,3152,4928,5664,6304};
    int iagpx[] = {0,640,1008,1600,2192,2560,3200};
    int iseq, ix1, ix2, ipx, irng, step[] = { 1,2,3,3,2,1};
    char ix, attrname[50];
    double sind = 0.0003104;
    float rad_2_deg = 180. / acos( -1 );
    int get_wt( float, float *, int, float *, char * );

    if (mside < 0) mside = 0;
    if (mside > 1) mside = 1;

    for (ib=0; ib<nbands; ib++) {
//        l1rec->polcor[ip*NBANDS+ib] = 1.0;
        l1rec->polcor[ip*nbands+ib] = 1.0;
    }

    if (l1rec->input->pol_opt == 0 || l1rec->input->pol_opt > 99)
        return;

    if (firstCall) {

      char  name   [H4_MAX_NC_NAME]  = "";
      char  sdsname[H4_MAX_NC_NAME]  = "";
      char  file   [FILENAME_MAX] = "";
      int32 sd_id;
      int32 sds_id; 
      int32 rank; 
      int32 nt; 
      int32 dims[H4_MAX_VAR_DIMS]; 
      int32 nattrs;
      int32 start[3] = {0,0,0}; 
      int32 end  [3] = {1,1,1}; 
      int32 status;

      int32  ndets;
      int32  nmside;
      int32_t  im, ia, id, ixb;


      for (ib=0; ib<nbands; ib++) {
          have_xcal[ib] = 0;
      }
      for (ixb=0; ixb<l1rec->input->xcal_nwave; ixb++) {
          if ((l1rec->input->xcal_opt[ixb] & XCALPOL) != 0) {
              //20180823
//	      if ((ib = bindex_get(l1rec->input->xcal_wave[ixb])) < 0) {
//	          HY1C_out("-E- %sline %d: xcal wavelength %f does not match sensor\n",
//	   	       __FILE__,__LINE__,l1rec->input->xcal_wave[ixb]);
//	          exit(1);
//	      };
              have_xcal[ib] = 1;
	  }
      }


      /* load polfile for each sensor wavelength */

      int wave[16];
      wave[0]=412;
      wave[1]=443;
      wave[2]=469;
      wave[3]=488;
      wave[4]=531;
      wave[5]=547;
      wave[6]=555;
      wave[7]=645;

      wave[8]=667;
      wave[9]=678;
      wave[10]=748;
      wave[11]=859;
      wave[12]=869;
      wave[13]=1240;
      wave[14]=1640;
      wave[15]=2130;
      for (ib=0; ib<nbands; ib++) {

//          sprintf(file,"%s%s%d%s",l1rec->input->polfile,"_",l1rec->iwave[ib],".hdf");
          sprintf(file,"%s%s%d%s",l1rec->input->polfile,"_",wave[ib],".hdf");
//   20180824
//          HY1C_out("Loading polarization file %s\n",file);

//          /* Open the file */
//          sd_id = SDstart(file, DFACC_RDONLY);
//          if (sd_id == FAIL){
//              fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
//              __FILE__,__LINE__,file,DFACC_RDONLY);
//              exit(1);
//          }

//          status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Detectors"),&ndets);
//          if (status != 0) {
//              HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
//                  __FILE__,__LINE__,"Number of Detectors",file);
//              exit(1);
//      }

//          status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Mirror Sides"),&nmside);
//          if (status != 0) {
//              HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
//                  __FILE__,__LINE__,"Number of Mirror Sides",file);
//              exit(1);
//      }
	  
//          if (ib == 0) {

//              // decide on which angle information needs to be read

//              if( sensorID == VIIRS ) {
//                strcpy( attrname, "Number of Scan angles" );
//                strcpy( sdsname, "scanangle" );
//              } else {
//                strcpy( attrname, "Number of AOIs" );
//                strcpy( sdsname, "AOI" );
//              }
//              status = SDreadattr(sd_id,SDfindattr(sd_id,attrname),&nang);
//              if (status != 0) {
//                  HY1C_out(
//                "-E- %s Line %d:  Error reading global attribute %s from %s.\n",
//                      __FILE__,__LINE__, attrname, file);
//                  exit(1);
//          }

//              if ((ang = (float *) malloc(nang*sizeof(float))) == NULL) {
//                  HY1C_out("-E- %s Line %d:  Error allocating memory for ang.\n",
//                      __FILE__,__LINE__);
//                  exit(1);
//          }
 
//              sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
//              status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
//              status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) ang);
//              if (status != 0) {
//                  HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
//                      __FILE__,__LINE__,sdsname,file);
//                  exit(1);
//              } else {
//                  status = SDendaccess(sds_id);
//              }
//      }


//          strcpy(sdsname,"m12");
//          sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
//          for (im=0; im<nmside; im++) for (ia=0; ia<nang; ia++) for (id=0; id<ndets; id++) {
//              start[0]=im; start[1]=ia; start[2]=id;
//              status = SDreaddata(sds_id, start, NULL, end, (VOIDP) &m12[ib][im][ia][id]);
//              if (status != 0) {
//                  HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
//                      __FILE__,__LINE__,sdsname,file);
//                  exit(1);
//              }
//          }
//          status = SDendaccess(sds_id);

//          strcpy(sdsname,"m13");
//          sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
//          for (im=0; im<nmside; im++) for (ia=0; ia<nang; ia++) for (id=0; id<ndets; id++) {
//              start[0]=im; start[1]=ia; start[2]=id;
//              status = SDreaddata(sds_id, start, NULL, end, (VOIDP) &m13[ib][im][ia][id]);
//              if (status != 0) {
//                  HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
//                      __FILE__,__LINE__,sdsname,file);
//                  exit(1);
//              }
          }
//          status = SDendaccess(sds_id);
//          status = SDend(sd_id);

          /* we need to translate from actual detectors to replicated detectors or  */
          /* aggregated detectors, to handle the HMODIS resolution switching. Get   */
          /* the conversion factor by comparing sensor ndets to table ndets.        */
	  if( sensorID == VIIRS )
            detfac[ib] = 1.;
          else
            detfac[ib] = l1rec->ndets * 1.0/ndets; /* 0.25 to 4 */
      }
      

//      /* get the max # pixels */
//      if (sensorID == HMODISA || sensorID == HMODIST) {
//	  switch (l1rec->input->resolution) {
//          case 250:
//  	    maxpix = 5416.0;
//            break;
//  	  case 500:
//  	    maxpix = 2708.0;
//            break;
//  	  default:
//  	    maxpix = 1354.0;
//            break;
//  	}
//      } else if (sensorID == MODISA || sensorID == MODIST) {
//          maxpix = 1354.0;
//      } else if (sensorID == VIIRS ) {
//          margin_s = l1rec->margin_s;
//          if( l1rec->scn_fmt == 0 )
//            maxpix = 3200.;
//          else
//            maxpix = 6304 + margin_s * 2;
//      } else {
// ////          HY1C_out("-E- %s line %d: Mirror geometry unknown for %s\n",
// ////                 __FILE__,__LINE__,sensorDir[sensorID]);
// //         exit(1);
//      }

      maxpix=2048;

      // precompute the weights and angle (AOI on mirror or scan angle)
      // dimension index for each pixel in the scan

      if( ( pxwt = (float *) malloc( maxpix * sizeof( float ) ) ) == NULL ) {
        HY1C_out("-E- %s Line %d:  Error allocating memory for weights.\n",
          __FILE__,__LINE__);
//        exit(1);
      }
      if( ( pxangix = (char *) malloc( maxpix * sizeof( char ) ) ) == NULL ) {
        HY1C_out("-E- %s Line %d:  Error allocating memory for angle index.\n",
          __FILE__,__LINE__);
//        exit(1);
      }

      // depending on instrument, get the angle and find the weight and index
      // into pol table for all pixels
ix=0;
//      if( sensorID == VIIRS ) {
//        if( l1rec->scn_fmt == 0 )   /* aggregated */
//          {
//          for( irng = 0; irng < 6; irng++ )
//            {
//              ix=ix;
////            ix1 = irng;
////            ix2 = irng + 1;
////            for( ipx = iagpx[ix1]; ipx < iagpx[ix2]; ipx++ ) {
////              iseq = ipx - iagpx[ix1];
////              zang = iagsm[ix1] + ( ( 2 * iseq + 1 ) * step[irng] - 1. ) / 2.;
////              zang = ( zang - 3151.5 ) * sind * rad_2_deg;
  
////              // use zang in same way as below - I hate to repeat this though
////              get_wt( zang, ang, nang, &wt, &ix );
////              pxwt[ipx] = wt;
////              pxangix[ipx] = ix;
////              }
//            }
//          }
//        else    /* unaggregated */
//          {
//          for( ipx = 0; ipx < maxpix; ipx++ )
//            {
////            zang = ( (float) ipx - 3151.5 - margin_s ) * sind * rad_2_deg;
////            get_wt( zang, ang, nang, &wt, &ix );
////            pxwt[ipx] = wt;
////            pxangix[ipx] = ix;
//              ipx=ipx;
//            }
//          }

//      else {

//        for( ipx = 0; ipx < maxpix; ipx++ ) {
////          zang = 10.5 + ipx / maxpix * (65.5 - 10.5);
////          get_wt( zang, ang, nang, &wt, &ix );
////          pxwt[ipx] = wt;
////          pxangix[ipx] = ix;
//                ix=ix;
//        }
//      }
//}
      nir_s  = bindex_get(l1rec->input->aer_wave_short);
      nir_l  = bindex_get(l1rec->input->aer_wave_long );
      
      firstCall = 0;


    /* apply sensitivites to radiances to get corrections */
//  20180824
    for (ib=0; ib<nbands; ib++) {

        ipb   = ip*nbands+ib;
//        idet  = (int32_t) rint(detnum/detfac[ib]);
//        alpha = l1rec->alpha[ip]/RADEG;
//        L_x   = l1rec->Lt[ipb]/l1rec->tg_sol[ipb]/l1rec->tg_sen[ipb];

//        if (L_x > 0.0) {
//	    L_qp = l1rec->L_q[ipb]*cos(2*alpha) + l1rec->L_u[ipb]*sin(2*alpha);
//	    L_up = l1rec->L_u[ipb]*cos(2*alpha) - l1rec->L_q[ipb]*sin(2*alpha);
// 	    if (have_xcal[ib]) {
//	        dm12 = get_xcal(l1rec,XM12,l1rec->iwave[ib]);
//	        dm13 = get_xcal(l1rec,XM13,l1rec->iwave[ib]);
//	        polcor[ipb] = 1.0 / ( 1.0 - dm12[ip] * L_qp/L_x - dm13[ip] * L_up/L_x );
//	    } else {
//	        for (iang=pxangix[pixnum]; iang<=(pxangix[pixnum] + 1); iang++) {
// 	            m1[iang] = 1.0 / ( 1.0
//                         - m12[ib][mside][iang][idet] * L_qp/L_x
//	                 - m13[ib][mside][iang][idet] * L_up/L_x );
//                }
//	        polcor[ipb] = m1[(int)pxangix[pixnum]] * (1.-pxwt[pixnum]) + m1[(int)pxangix[pixnum] + 1] * pxwt[pixnum];
//	    }
//            dpol  [ipb] = sqrt(pow(l1rec->L_q[ipb],2.0)+pow(l1rec->L_u[ipb],2.0))/L_x;

            /* quick-fix to polcor of aerosol bands only */
        if (l1rec->input->pol_opt == 6 && ib != nir_l && ib != nir_s){
	        polcor[ipb] = 1.0;
            }
//        } else {
	    polcor[ipb] = 1.0;
         dpol  [ipb] = 0.0;
//	}
//    }
}
}


int get_wt( float zang, float *ang, int nang, float *wt, char *ix )
/*******************************************************************

   get_wt

   purpose: derive a weight and index into angle array

   Returns type: int - 0 if all is OK

   Parameters: (in calling order)
      Type              Name            I/O     Description
      ----              ----            ---     -----------
      float             zang             I      angle to interpolate to
      float *           ang              I      array of angles in 
                                                interpolation table
      int               nang             I      # of angle tie points
      float *           wt               O      derived weight
      char *            ix               O      index of first angle to use in 
                                                weighting

*******************************************************************/
  {
  int iang, ix2;
  
  if (zang <= ang[0])
    *ix = 0;
  else if (zang >= ang[nang-1])
    *ix = nang - 2;
  else
    {
    for (iang=1; iang<nang; iang++)
      {
      if (zang < ang[iang])
        {
        *ix = iang - 1;
        break;
        }
      }
    }
  ix2 = *ix + 1;
  *wt = ( zang - ang[(int)*ix] ) / ( ang[ix2] - ang[(int)*ix] );
  return 0;
  }
